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Abstract 

When a block made of an elastomer is subjected to large shear, 
its surface remains flat. When a block of biological soft tissue is sub- 
jected to large shear, it is likely that its surface in the plane of shear 
will buckle (apparition of wrinkles). One factor that distinguishes 
soft tissues from rubber-like solids is the presence - sometimes visi- 
ble to the naked eye - of oriented collagen fibre bundles, which are 
stiffer than the elastin matrix into which they are embedded but are 
nonetheless flexible and extensible. Here we show that the simplest 
model of isotropic nonlinear elasticity, namely the incompressible neo- 
Hookean model, suffers surface instability in shear only at tremendous 
amounts of shear, i.e., above 3.09, which corresponds to a 72° angle 
of shear. Next we incorporate a family of parallel fibres in the model 
and show that the resulting solid can be either reinforced or strongly 
weakened with respect to surface instability, depending on the angle 
between the fibres and the direction of shear, and depending on the 
ratio E/jj, between the stiffness of the fibres and that of the matrix. 
For this ratio we use values compatible with experimental data on soft 
tissues. Broadly speaking, we find that the surface becomes rapidly 
unstable when the shear takes place "against" the fibres, and that as 
E/fj, increases, so does the sector of angles where early instability is 
expected to occur. 
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1 Introduction 

Rubber-like solids and biological soft tissues can both be efficiently mod- 
elled within the framework of finite elasticity, which can account for large 
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deformations, physical nonlinearities, incompressibility, residual stresses, vis- 
coelasticity, etc. One of the most salient differences between the two types 
of solids is that at rest, elastomers are essentially isotropic whilst soft tissues 
are essentially anisotropic, because of the presence of collagen fibre bundles. 
In that respect, it is worthwhile to consider the effect of incorporating fam- 
ilies of parallel fibres into an isotropic matrix, and see if it can model some 
striking differences between the mechanical behaviour of elastomers and of 
soft tissues. Consider for instance the large shear of a solid block. When 
the block is made of an elastomer such as silicone, its surface remains sta- 
ble; when it is made of a biological soft tissue such as skeletal muscle, its 
surface wrinkles for certain ranges of orientation between the direction of 
shear and the (presumed) direction of fibres, see Fig. [TJ Here we show that 
one of the simplest models of anisotropic nonlinear elasticity, which requires 
only knowledge of the fibre/matrix stiffness ratio, is sufficient to successfully 
predict these behaviours. 

To model the isotropic elastomer (Section 2), we take the incompressible 
neo-Hookean solid, and find that it does not suffer surface instability unless it 
is subjected to a substantial amount of shear (critical amount of shear: 3.09, 
critical angle of shear: 72°). In that case the wrinkles are aligned with the 
direction of greatest stretch. (The wrinkling analysis relies on the incremental 
theory of nonlinear elasticity, see for instance Biot [1] or Ogden [2]). Next, 
we introduce one family of parallel fibres into the model (Section 3). To 
model biological soft tissues with one preferred direction (Section 4), we take 
the incompressible neo-Hookean strain energy density, augmented by the so- 
called 'standard reinforcing model': this model has only two parameters, 
namely the shear modulus /i of the soft (neo-Hookean) matrix and the fibre 
stiffness E. 

With respect to surface instability, only the ratio E/fi of these two quan- 
tities plays a role. We take it to be equal in turn to 40.0, 20.0, and 10.0, 
in agreement with the range of experimental measures found in the litera- 
ture. We then find that when the angle between the direction of shear and 
the direction of the fibres is small, the solid is much more stable than the 
isotropic solid obtained in the absence of fibres; when the angle increases 
but is less than 99.0° (for E/fi = 40.0), 102.8° (for E/fi = 20.0), 108.1° (for 
Ej fi = 10.0), the solid remains more stable than the isotropic solid; however, 
when the angle exceeds those values, the critical amount of shear for surface 
instability drops to extremely low levels, indicating the appearance of wrin- 
kles as soon as shearing occurs. In that case, the wrinkles are found to be 
almost orthogonal to the fibres, in accordance with visual observations. 

It is hoped that the paper provides a greater understanding of the causes 
of certain instabilities in soft tissues and a quantitative tool to measure what 
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Figure 1: Shearing (along the arrows) a block of silicone of approximate size 
15cm x 10cm x 1.5cm and a block of mammalian skeletal muscle (beef) of 
approximate size 15cm x 10cm x 3cm; one does not exhibit surface instability, 
the other does. 

deformations (critical amounts of shear) are permissible and in which direc- 
tions. Surface instability has a direct connection to slab and tube buckling, 
which in biomechanics may potentially translate into aneurysms formation, 
arterial kinking and tortuosity, brain trauma, and many other, still not well 
understood, pathologies. 

2 Surface instability of a sheared isotropic 
solid 

First, we recall known results in the theory of surface wrinkling valid for 
isotropic solids. 

Consider a semi-infinite body made of an incompressible isotropic neo- 
Hookean solid, for which the strain energy function W, written as a function 
of the principal stretch ratios Ax, A 2 , A 3 , is given by 

W = /i (A? + A! + Al-3)/2. (1) 

Here /i is the shear modulus, and A1A2A3 = 1 by the incompressibility con- 
straint. Then subject the solid to a large homogeneous static deformation, 
such that A2 is the stretch ratio along the normal to the free surface. It 
has long been known that the surface becomes unstable when the following 

wrinkling condition is met, 

A?A 3 = o-q, (2) 

where <To — 0.296 is the real root of a 3 + a 2 + 3a — 1 = (Green and Zerna [3], 
Biot [1]). 
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Figure 2: Large plane strain deformation of a unit cube near the surface of a 
semi-infinite incompressible neo-Hookean solid. When the solid is compressed 
by 71% (or equivalently, stretched by 238%), its surface wrinkles. Note that 
the analysis quantifies neither the amplitude nor wavelength of the wrinkles. 

In the following plane strain situation, 

Ai = A, A 2 = l, A 3 = A-\ (3) 

the critical stretch of compression found from Eq. is clearly Ai = cr — 
0.296 (and then A3 = (Tq 1 ~ 3.38). The conclusion is that when a semi-infinite 
neo-Hookean solid, which is neither allowed to expand nor contract along 
the normal to its boundary, is compressed by 71% in a given direction (lying 
in the boundary), it buckles with wrinkles developing along the direction 
orthogonal to the direction of compression. Equivalently, when it is stretched 
by 238%, it buckles with wrinkles parallel to the direction of tension. Figured 
summarizes these results. 

It is natural to wonder whether the surface might have become unstable 
in other directions earlier, that is at compressive (^ 1) ratios larger than 
0.296, or at tensile (^ 1) ratios smaller than 3.38. Flavin [4] shows that 
wrinkles develop parallel to the direction making an angle 9 with the principal 
direction of strain associated with the stretch ratio A3 when the following 
wrinkling condition is met 

A 2 A 2 (A 2 cos 2 # + A 2 sin 2 #) = a 2 . (4) 

In the plane strain situation Eq. ([3]), this condition is quadratic in A 2 , 

A 4 cos 2 0- A 2 t^ + sin 2 = O. (5) 
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Figure 3: Large simple shear of a unit square in the surface of a semi-infinite 
incompressible neo-Hookean solid. When the solid is sheared by an amount 
Kq ~ 3.09 (Figure on the right), its surface wrinkles. The corresponding 
angle of shear is tan -1 K Q ~ 72.0°, which physically, is abnormaly large. 
Then the wrinkles are parallel to the direction of greatest tension, which 
makes an angle (po — 16.5° with the direction of shear (and so, the wrinkles 
are almost aligned with the sheared faces.) 

It has real roots provided 9 is in the ranges —9q ^ 9 ^ 9$ oi n /2 — 9q ^ 
9 ^ 7r/2 + 9q, where 9$ = (1/2) sin -1 o"q ~ 2.51°. In the former range, the 
compressive critical stretch found from the biquadratic Eq. §5§ turns out to 
be smaller than 0.296 and in the latter range, to be larger than 3.38. Thus 
surface instability for plane strain Eq. ([3]) occurs when the isotropic neo- 
Hookean half-space is in compression at a ratio o~o or equivalently, in tension 
at a ratio er^ 1 . The wrinkles are parallel to the direction of greatest stretch 
and orthogonal to the direction of greatest compression. 

Now simple shear belongs to the family of plane strains Eq. ([H]), with the 
following connection between the principal stretches and the amount of shear 
K (see Ogden [2]) for instance), 

K = \-\~\ A = K/2 + a/1 + K 2 /4. (6) 

Also, the direction of greatest stretch is at an angle ip with the direction of 
shear, where ip G]0,7r/4] is given by 

tan2^ = 2/AT. (7) 

Clearly A > 1 here, and so surface shear instability occurs in tension, when 
the amount of shear is equal to K = a x — o~o — 3.09. The corresponding 
critical angle of shear is then tan -1 Kq ~ 72.0°, see Fig. [31 This is quite large 
shear. 
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3 Sheared fibre-reinforced solids 



3.1 Finite simple shear 

Now we consider a semi-infinite composite incompressible solid, made of an 
isotropic matrix reinforced with one family of parallel extensible fibres, them- 
selves parallel to the boundary of the solid. In the undeformed configuration, 
we call (Xi,X 2 ,X 3 ) the set of Cartesian coordinates such that the solid is 
located in the X 2 ^ region. We denote by Ei, E 2 , E 3 the orthogonal unit 
vectors defining the Lagrangian (reference) axes, aligned with the X\, X 2 , 
X 3 directions, respectively. 

When the solid is sheared in the direction of E±, the particle at X moves 
to its current position x. We call F = dx/dX the associated deformation 
gradient tensor, and B = FF T the left Cauchy-Green strain tensor. We 
then call (xi,x 2 ,x 3 ) the Cartesian coordinates, aligned with (X 1? X 2 , X 3 ), 
corresponding to the current position x. In the current configuration, the 
basis vectors are e±, e 2 , e 3 , and here they are such that e,; = Ei {i = 1, 2, 3). 
The simple shear of amount K is described by 

x 1 = X 1 + KX 3 , x 2 = X 2 , x 3 = X 3 . (8) 

We thus find in turn that 

F = I + Ke 1 <g> E 3 , B = I + K{e l ®e 3 + e 1 ®e 3 )+K 2 e 1 ®e l . (9) 

The principal stretches are given by Eq. ([3]) and Eq. (Q, and the first principal 
isotropic invariant l\ = tr B is given here by 

J x = 3 + K 2 . (10) 

Note that for shear, the second principal isotropic invariant, I 2 = [If — 
tr (B 2 )}/2 is also equal to 3 + K 2 . 

3.2 One family of fibres 

For solids reinforced with one family of parallel fibres lying in the plane of 
shear, we work in all generality and consider that the angle $ (say) between 
the fibres and the X± direction can take any value. In other words, the unit 
vector M (say) in the preferred fibre direction has components 

M = cos$£i + sin $£ 3 , (11) 
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Figure 4: A unit square lying on the surface of a semi- infinite solid reinforced 
with one family of fibres (thin lines) and subject to a simple shear of amount 
K = 0.5 (angle of shear is tan -1 K ~ 26.6°) in the X± direction. In the 
reference configuration, the fibres are along the unit vector M, at the angle 
$ = 60° with the Xx-axis. In the current configuration, they are along m. 
The unit vector n is orthogonal to the wrinkles' front (when they exist). 
Finally, the dashed line is aligned with the direction of greatest stretch; it is 
at an angle ip ~ 38° to the direction of shear. 



in the reference configuration. Simple shear is a homogeneous deformation, 
and so M is transformed into m = FM in the current configuration, that 
is 

m = (cos$ + Ksm$)e 1 + sin$e 3 . (12) 

Without loss of generality, we take the ranges K ^ 0, ^ $ ^ n, which 
cover all possible orientations of the fibres with respect to the direction of 
shear. 

To fix the ideas, consider Fig. HI There we shear the half-space by a 
finite amount K = 0.5, in the direction making an angle $ = 60° with the 
fibres. Notice that a unit vector n making an angle 9 with the direction of 
shear is also represented in the current configuration. This is the normal to 
the wrinkles' front; in the next section we look for surface wrinkles in all 
directions (the angle 9 spans the interval [0°, 180°]) and we determine which 
is the smallest corresponding critical amount of shear. 

Finally we introduce the anisotropic invariants I4 = m • m and I5 = 
Fm-Fm; in particular we find 

h = 1 + Ksm2$ + K 2 sin 2 $. (13) 

Recall that I4 is the squared stretch in the fibre direction [5]. In particular, 
if I4 ^ 1 then the fibres are in extension, and if I4 ^ 1 then they are in 
compression. Clearly here, when ^ $ ^ 7r/2, the fibres are always in 



7 



extension but when 71/ 2 < $ < tt, there exist a certain amount of shear 
(explicitly, —2/ tan$) below which the fibres are in compression. 



3.3 Constitutive assumptions 

In general, the strain-energy density W of a hyperelastic incompressible 
solid reinforced with one family of parallel extensible fibres depends on the 
isotropic invariants I\ and I2, and on the anisotropic invariants (Spencer [5]) 
1^ and I5. We assume that W is the sum of an isotropic part and an 
anisotropic part. For the isotropic part, modelling the properties of the 'soft' 
matrix, we take the neo-Hookean strain-energy density in order to make a 
connection with the results of Section 2. For the anisotropic part, modelling 
the properties of the extensible 'stiff' fibres, we take a function of I± only, 
say F(I±). Hence, we restrict our attention to those solids with strain energy 
density 



This assumption is quite common in the biomechanics literature. Although 
it does not prove crucial to the analysis, it leads to compact and revealing 
expressions (Note that the consideration of a more general W poses no major 
extra difficulty, but results in much longer expressions.). 

The corresponding Cauchy stress tensor <x is (see e.g. [6]): a = —pi + 
p,B + 2F'(J 4 )m ® m, where p is a Lagrange multiplier introduced by the 
constraint of incompressibility. The surface x 2 — is free of tractions: here 
012 = °23 = follows from B\ 2 = B23 = and m • e 2 = (see Eq. ([9]) 
and Eq. (fi"2j) ). whilst cr 2 2 = gives p — fi. Thus, the pre-stress necessary to 
maintain the shear Eq. (jSJ) is 



showing that the directions of principal stress and strain do not coincide in 
general (except when the preferred direction is aligned with principal direc- 
tions of strain). 

4 Surface instability 

4.1 Incremental deformations 

We seek solutions to the incremental equations of equilibrium and incre- 
mental boundary conditions in the form of a sinusoidal perturbations whose 
amplitude decays rapidly with depth. In contrast to the isotropic case of 




(14) 




(15) 



8 



Section 2, we do not know a priori in which direction the wrinkles should be 
aligned, and we take the normal to the wrinkles' front n (say) to lie in the 
(xxXs) plane at an arbitrary angle 9 with x±, see Fig. HJ Hence, we seek a 
perturbation solution u (mechanical displacement) and p (increment of the 
Lagrange multiplier associated with incompressibility) in the form, 



ife(cos 0zi+sin 8x3) 



{u,p} = {U(kx 2 ),ikP(kx 2 )}e 



where k is the "wave" -number and U, P are functions of kx 2 alone. 
The incremental equations read 



(16) 



0. 



u 



.1 -j 



0. 



(17) 



where the comma denotes partial differentiation with respect to Xj, and s is 
the incremental nominal stress tensor. Its components are [2], 



= A jUkUk,i + pu j:i - p8, 



where *4o is the fourth-order tensor of instantaneous elastic moduli. In gen- 
eral it has a long expression for fibre-reinforced solids, with possibly 45 non- 
zero components, see for example [7,8]. For W in the form Eq. ( |T4l) . B by 
Eq. and M by Eq. (|TT|) . we find the following components 



A 



Ojilk 



^5 ik Bji + 2F'(I 4 )5 ik m j mi + 4F"( y I 4 )m i m j m k mi 



(19) 



see Merodio and Ogden [9]. Clearly, these components have the symmetries 
Aojiik = Aoikji and Aojuk = Aojkii- We end up with 23 non-zero components, 
several of which are equal to one another (in toto there are 13 different 
components). 

Clearly, if u and p are of the form Eq. ( TL6l) . then by Eq. ( TL8l) the Sji are 
of a similar form, say 



ikSji(kx 2 )e 



ife(cos Srri+sin 6x3) 



(20) 



where the Sji are functions of the variable kx 2 only. By a systematic proce- 
dure, first laid down by Chadwick [10] (see also [11-14]), we can eliminate P 
and write the incremental equations of equilibrium as a first-order differential 
system. This is known as the Stroh formulation of the problem, 



U' 
S' 



iN 



U 

s 









S21 


where U = 


u 2 


, s = 


s 22 




u 3 




5*23 



N 



N x N 2 
N 3 N 1 



(21) 
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and the symmetric 3x3 matrices N±, N 2 , N 3 are given by 

rj k 
v 
k x 
(22) 
with 

77 = (-4.011H + 3/i) cos 2 + 2^4o 1131 cos 6 1 sin + 4 3i3i sin 2 9, 

v = 4 m2 cos 2 9 + 2^01232 cos 9 sin # + A)3232 sin 2 # - /i, 

X = A1313 cos 2 6» + 2^01333 cos 9 sin 6> + (A3333 + 3/i) sin 2 9, 

K = Ami3 cos 2 9 + (2.4.01133 + 3/i) cos 9 sin # + A)3i33 sin 2 6 1 . (23) 

Notice how all the information relative to anisotropy is located in the iV 3 
matrix. 

The solution to the system Eq. f[2"Tj) is clearly an exponential 

{U,S} = {U°,S°}e ikqX2 , (24) 

where U°, S° are constant vectors and q is an eigenvalue of N. The charac- 
teristic equation associated with N is a bicubic [8], 

g 6 -(2-^V + (i + ^ + ^V + i^ = o, ( 25 ) 
V n J V / n 

where the quantity e is defined by 

e = x cos 2 9 — 2k cos # sin # + 77 sin 2 9. (26) 

The existence of real roots to this equation corresponds to the loss of ellip- 
ticity of the governing equations (material instabilities). This possibility has 
been thoroughly investigated before, see [9,15,16]. Here we focus on complex 
roots and keep those satisfying Img > , for a surface-type bifurcation which 
decays with depth (geometric instability). 

4.2 Wrinkling condition and resolution scheme 

Over the years, many schemes have been developed to solve surface boundary 
problems using the Stroh formulation; we used in turn the determinantal 
method [17], the Riccati matrix equation of surface impedance [13,14], and 
explicit polynomial equations [18], in order to double-check our numerical 
computations. 
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cos 9 
-iVi = cos6» sin# 
sin^ 



No 



1/fi 

l/fi 



-N, 



The crucial boundary condition is to find the amount of shear at which 
the surface of the sheared solid is free of tractions. The safest way to express 
this is 

det Z = 0, (27) 

where Z is the (Hermitian) surface impedance matrix, which relates tractions 
to displacements through S = iZU. We remark that the schemes are not as 
safe in surface stability problems as they are in surface wave theory because 
of incompressibility [13, 14] and non-monotonicity of det Z with K. 

Once Eq. (|27[) is reached, we can construct an incremental solution to the 
equations of equilibrium which is adjacent to the large shear equilibrium, and 
signals the onset of surface instability. We adopted the following strategy: 

(i) Fix $, the angle between the direction of shear and the preferred di- 
rection; 

(ii) Fix 9, the angle between the direction of shear and the normal to the 
wrinkles' front; 

(Hi) Find (if it exists) the corresponding critical amount of shear such that 
Eq. f[2"Tj) is satisfied. 

Then repeat Steps (ii) and (Hi) for other angles 9 until the entire surface 
is spanned, and keep the smallest critical amount of shear K cr (say) for the 
angle $ chosen in Step (i). Then take a different value of $, until all possible 
fibre orientations are covered. In fine a graph of K cr as a function of $ is 
generated. 

5 Numerical results for biological soft tissues 

We take the standard reinforcing model, 

W = n(I x - 3)/2 + E(h - l) 2 /4, (28) 

where E is an extensional modulus in the fibre direction. This model has been 
used for several soft tissues, such as papillary muscle [19], myocardium [19], 
skeletal muscles [20], or brainstem [21]. 

That latter reference examines the ability of the constitutive model Eq. ( |28l) 
to describe the mechanical response of porcine brainstem specimens. Recall 
that large deformations, in particular large shears, of brain tissue are often 
associated with traumatic brain injuries (Doorly and Gilchrist, 2006). Ning 
et al. [21] find that the model provides good agreement with experimental 
data; they estimate that for 4 week old pigs, E is about 20 times larger than 
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Figure 5: Variations of the critical amount of shear for surface instability 
with the angle between the directions of shear and the fibres. The solid 
is modelled as a neo-Hookean matrix reinforced with one family of fibres 
(standard reinforcing model); the ratio of the matrix shear modulus to the 
fibre stiffness is taken in turn as 40.0, 20.0, and 10.0. The 3 graphs coincide 
as long as < $ < $ , where $ = 99.0°, 102.8°, 108.1°, respectively. At 
$ ~ $0) the half-space switches from being very stable (K cr > 3.09) to being 
easily unstable (K cr < 0.3). The part of the plot corresponding to K CT > 5 
is not shown for physical and visual reasons. 



fj,. In a recent review on physical properties of tissues for arterial ultrasound, 
Hoskins [23] emphasizes the need for constitutive models of nonlinear elastic 
behavior. He also collects available data for arterial walls: in particular for 
abdominal aortic aneurysms, ex vivo measurements indicate that E is about 
10 times larger than // whilst for human atherosclerotic plaque, E seems to 
be more than 40 times /i. For our numerical computations we take in turn 
the values E/fi = 40.0, 20.0, 10.0, and collect the corresponding results on 
Fig. E 

Broadly speaking, we find a region where the solid is strongly reinforced 
by the family of fibres, followed by an abrupt drop in the value of the critical 
amount of shear for surface instability, which occurs earlier as E//J, increases. 

When the fibres are aligned with the direction of shear, they are not 
stretched and they play no role; thus it is appropriate that at $ = 0.0°, we 
find K cr = 3.09, the critical amount of shear for an isotropic neo-Hookean 
half-space, see Section 2. 

Next we find that K cr shoots up to unrealistic values when <£> > 0.0°: for 
instance K CT = 32.48 when $ = 3.0° (not represented for visual convenience). 
Hence, the solid is strongly reinforced with respect to surface stability when 
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the shear takes place more or less along the fibres: wrinkling is prevented. 

As the angle $ between the shear and the fibres increases, the critical 
amount of shear goes through a maximum, then a minimum, always re- 
maining above 3.09, the value for an isotropic neo-Hookean half-space, as 
long as $ ^ $ , where $ = 99.0°, 102.8°, 108.1°, approximatively, for 
E/fj, = 40.0,20.0,10.0, respectively. It is worth noting that in the range 
90.0° < $ < $ , the fibres undergo a slight compression at low shear levels, 
and then are in extension until the critical amount of shear is reached; even 
when the fibres are compressed, the half-space remains stable. 

When the angle $ is large, $o < 3? < 180.0°, the half-space becomes 
unstable at low amounts of shear. For instance at $ = 99.07°, we find that 
K CT = 0.153 when E/fi = 40.0; note that in reaching that critical amount 
of shear, the fibres are compressed by less than 1.3%. The switch from high 
to low critical amounts of shear is abrupt, due to the non-monotonicity of 
det Z with K: this quantity has a minimum in the high range [K > 3.09) 
which is always negative (indicating the existence of a root to Eq. (|27|) ). but 
it can also have a minimum in the low range (K < 0.3). This minimum is 
positive when $ < <£> ( no r °°t to Eq. (12T1) ) but negative when $ > $ , hence 
the jump in K CT . 

Finally we note that in the range <3?o < $ < 180.0°, the angle 9 normal 
to the wrinkles' front is close to $ (within 2°), indicating that the wrinkles 
are almost at right-angle with the fibres; these predictions are in accordance 
with the observation of Fig. [TJ 

6 Discussion 

We developed a quantitative methodology to understand the formation of 
wrinkles in some biological soft tissues. The analysis allowed us to model 
some visual observations of a sheared elastomer versus a sheared piece of 
skeletal muscle, based on a simple nonlinear anisotropic constitutive law (re- 
quiring the knowledge of only one quantity, Ej \x). 

Studying the geometry and mechanics of wrinkles is relevant to many 
biomechanical applications such as for instance the treatment of scars, and 
our results may provide some help in developing rational approaches to these 
problems. The next logical step is to apply and generalize this methodology 
to model the wrinkling of skin and other biological membranes. These may 
require more work than here, with the consideration of two families of parallel 
fibres (the collagen network), but the methodology remains essentially the 
same. It is also exact, versatile, and more convenient to apply than methods 
based on approximate theories (e.g. Foppl-von Karman plate equations) 
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because it can accommodate easily anisotropy, nonlinear constitutive laws, 
finite thickness, and large homogeneous pre-deformation. 
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